Garrett ef al. BMC Medical Genomics 2013, 6(Suppl 1):S2 
http://www.biomedcentral.eom/1755-8794/6/S1/S2 



Medical Genomics 



RESEARCH Open Access 



Short and long term gene expression variation 
and networking in human proximal tubule cells 
when exposed to cadmium 

Scott H Garrett 1 , Kaitlin Clarke 2 , Donald A Sens 1 , Youping Deng 3 , Seema Somji 1 , Ke K Zhang 1,2 * 

From The 2011 International Conference on Bioinformatics and Computational Biology (BIOCOMP'1 1) 
Las Vegas, NV, USA. 18-21 July 2011 



Abstract 

Cadmium (Cd 2+ ) is a known nephrotoxin causing tubular necrosis during acute exposure and potentially 
contributing to renal failure in chronic long-term exposure. To investigate changes in global gene expression 
elicited by cadmium, an in-vitro exposure system was developed from cultures of human renal epithelial cells 
derived from cortical tissue obtained from nephrectomies. These cultures exhibit many of the qualities of proximal 
tubule cells. Using these cells, a study was performed to determine the cadmium-induced global gene expression 
changes after short-term (1 day, 9, 27, and 45 uM) and long-term cadmium exposure (13 days, 4.5, 9, and 27 uM). 
These studies revealed fundamental differences in the types of genes expressed during each of these time points. 
The obtained data was further analyzed using regression to identify cadmium toxicity responsive genes. Regression 
analysis showed 403 genes were induced and 522 genes were repressed by Cd 2+ within 1 day, and 366 and 517 
genes were induced and repressed, respectively, after 13 days. We developed a gene set enrichment analysis 
method to identify the cadmium induced pathways that are unique in comparison to traditional approaches. The 
perturbation of global gene expression by various Cd 2+ concentrations and multiple time points enabled us to 
study the transcriptional dynamics and gene interaction using a mutual information-based network model. The 
most prominent network module consisted of INHBA, KIF20A, DNAJA4, AKAP12, ZFAND2A, AKR1B10, SCL7A11, and 
AKR1C1. 



Introduction 

Cadmium is a toxic heavy metal that is widely distributed 
due to industrial pollution such as mining, refining of 
metals, burning of fossile fuels, battery production and the 
use of tobacco. The health threats presented by cadmium 
has been increasingly recognized by government agencies, 
healthcare providers and the general public [1]. One of the 
primary target organs of cadmium toxicity is the human 
kidney where cadmium accumulates and induces a num- 
ber of adverse events. The proximal tubule is the main site 
for cadmium accumulation and the site of cadmium 
induced toxicity. When the levels of cadmium exceed a 
tolerance level of the intracellular defense systems, such as 
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metallothionein and glutathione, cells become susceptive 
to the toxic effect of this heavy metal. [2,3]. Under high 
dose acute exposure, cadmium is known to cause proximal 
tubule necrosis. On the other hand, a few cohort studies 
showed chronic low doses of cadmium can also induce 
proximal tubule toxicity [4,5]. Toxicokinetic models indi- 
cated that the half-life of cadmium in humans is between 
15 and 20 years [6]. This suggests that low level exposure 
of cadmium may be harmful to human kidney via the long 
term accumulation of this metal in proximal tubule. In 
this study, we investigate the acute and chronic effects of 
cadmium on human proximal tubule (HPT) cells. The 
HPT cell culture system developed in our lab has been 
shown to have important properties that resemble human 
proximal tubule [7] . The ability of forming dome by con- 
fluent cell monolayers implies the HPT cells retain the 
vectorial transport activity. The cells are mortal and 
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exposure to high concentration of cadmium induces a 
necrotic mechanism of cell death [8,9]. By exposing HPT 
cells to various doses of cadmium in either a short period 
of 1 day or a long period of 13 days, we expect different 
physiological responses to be elicited which represent the 
acute and chronic responses to cadmium toxicity. In order 
to investigate the underlying molecular mechanisms of cell 
toxicity and compensatory cellular defenses related to cad- 
mium exposure, we profiled the global gene expression in 
HPT cells exposed to increasing concentrations of this 
metal during short-term (1 day) and long-term (13 days) 
exposure. 

An initial analysis of the microarray data has previously 
reported a number of genes that were differentially 
expressed between cadmium treated and control samples 
[10]. In this manuscript, we investigated how genes 
respond to different doses of cadmium and how the gene 
network is formed in response to cadmium toxicity. The 
high cost of microarray chips prevents researchers from 
collecting microarray data with large sample sizes; thereby, 
limiting the power of statistical analysis. In this study, we 
only have one sample for each condition-time point, how- 
ever, the overall experimental design enabled us to investi- 
gate the gene expression profiles correlated with cadmium 
concentration. Furthermore, we developed a gene set 
enrichment analysis (GSEA) method that can be applied 
to regression-based tests for individual genes, while tradi- 
tional GSEA methods focus on two-sample comparison 
analysis, such as t test. 

This study adopted a complicated design by perturbing 
HPT cells with a variety of doses of cadmium and under 
various time periods that resulted in variation and 
dynamics of global gene expression. Therefore, it provides 
a unique opportunity to investigate the dynamic change of 
gene expression levels and the network structure of gene 
interactions in a well-designed experimental environment. 
A mutual information-based network model [11] was used 
to investigate the gene networking in response to cad- 
mium toxicity and a network module with a small set of 
activated genes with high interaction between each other 
was identified. 

Materials and methods 

Cell culture and total RNA collection 

The HPT cell culture system and cadmium treatment 
has been described previously [7,10]. RNA was purified 
from triplicate cultures of HPT cells by the RNeasy 
Mini kit (Qiagen, Valencia, CA). The integrity and qual- 
ity of RNA was ensured via the ribosomal bands using 
agarose gel electrophoresis. 

Microarray analysis 

Microarray experiment was performed by Genome 
Explorations Inc. (Memphis, TN) using Affymetrix 



GeneChip Human Genome U133 Plus 2.0 arrays, which 
cover over 47,000 human transcripts and variants repre- 
senting approximately 39,000 of the best characterized 
human gene transcripts and expressed sequence tags 
(ests) (Affymetrix, Santa Clara, CA). The scanned 
images were analyzed using programs resident in Gene- 
Chip Operating System vl.4 (GCOS; Affymetrix). The 
raw data was processed and normalized using MAS 5.0 
statistical algorithm that gave detection call regarding 
signal present or absent for each probe. The correlation 
of gene expression with cadmium concentrations was 
assessed using linear regression. The gene signal path- 
ways that were activated by cadmium were identified 
using gene set enrichment analysis (GSEA) [12]. The R 
package GSA that was based on standardized statistics 
was employed for gene set enrichment analysis [13]. 
The gene list for each pathway was obtained from the 
Kyoto Encyclopedia of Genes and Genomes (KEGG) 
http://www.genome.jp/kegg. Cluster analysis was used to 
group the microarray samples using hierarchical cluster- 
ing method. The distances between samples were calcu- 
lated by Pearson dissimilarity and the agglomeration of 
clusters was performed by Ward linkage method. 

Networks analysis 

The relations between genes were defined by mutual 
information [11]. The resulting adjacent matrix was used 
for networks analysis. Network modules were identified by 
hierarchical clustering. The set of genes that have short 
distances between each other and large distance with 
other genes was defined as a network module. The most 
prominent network module was identified for further 
investigation. All data analyses were performed using R 
statistical programming language. 

Results 

Gene expression variation in response to acute (1 day) 
cadmium toxicity 

The HPT cells were exposed to various levels of concen- 
trations of Cd 2+ for a short term of 24 hours. The Cd 2+ 
concentrations were 0, 9, 27, and 45 \a,M. Our hypothesis 
is that Cd 2+ exposure will elicit an acute toxicity response 
to HPT cells within one day. As reported before, the HPT 
cells exhibited a stress response during 24 hour exposure 
to increasing concentrations of Cd 2+ and showed visible 
signs of toxicity at the high dose [10]. The goal of the 
current study is to discover the genes that are responsive 
to the short-term exposure of Cd 2+ using messenger 
RNA microarray data. For each Cd 2+ concentration, a 
sample of HPT cells was harvested, total RNA was puri- 
fied and submitted for microarray analysis. The Affyme- 
trix U133 Plus 2.0 array was used for this analysis. When 
using MAS 5.0 for testing, we found 28,720 genes that 
were present in at least one sample. These genes were 
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used for subsequent analysis. In order to identify the 
genes whose expression levels were dependent on Cd 2+ 
concentrations, we conducted regression analysis for cor- 
relating gene expression with concentrations. This 
resulted in 4,036 genes that have significant P values (< 
0.05). However, none of them had a false discovery rate 
below 0.05 due to the small size of only 4 samples. 
Figure la shows the P value distribution that was tilted 
to low values implying that many genes had positive 
responses to acute Cd 2+ toxicity. If only a few genes 



had positive responses, the P values would be evenly 
distributed. Using a P value cutoff of 0.01, we reported 
403 genes that were positively correlated with cadmium 
concentration (Table 1), and 522 genes that were nega- 
tively correlated (Table 2). The positively correlated 
genes in the top list included heat shock genes like 
HSP40, HSP70 and HSPBAP1, and ion and anion trans- 
porters SLC7A11, SLC16A1, SLC20A2, SLC25A37 and 
SLC7A1. Both the highest positively and negatively cor- 
related genes contained cell cycle arrest and cell death 
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Figure 1 The distributions of P values in the regression analysis of gene expression versus cadmium concentration, (a) The histogram 
of P values for 1 day exposure of cadmium, (b) The histogram of P values for 13 days exposure of cadmium. 
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Table 1 Cd 2+ 


induced genes in acute cadmium-exposed human proximal tubule cells (1 day, p value < 0.001) 


Probe 


Symbol 


Gene name 


Pr(> |t|) 


22506 1_at 


DNAJA4 


DnaJ (Hsp40) homolog, subfamily A, member 4 


1.16E-05 


207528_s_at 


SLC7A1 1 


solute carrier family 7, (cationic amino acid transporter, y+ system) member 1 1 


1 .66E-05 


209774_x_at 


CXCL2 


chemokine (C-X-C motif) ligand 2 


4.83 E-05 


202744_at 


SLC20A2 


snlutp Crirripr familv 70 fnhosnh^tp trrinsnnrtpr.1 mpmhpr 7 


5 00F-05 


213524_s_at 


G0S2 


G0/G1 switch 2 


7.15E-05 


220533_at 


NA 


NA 


7.69E-05 


208544_at 


ADRA2B 


arl rpnprnir alnh^-^R- rprpntnr 


7.89E-05 


23921 1_at 


NA 


NA 


8.93E-05 


225339_at 


SPAG9 


enprm a^soristpH ^ntiripn Q 


0 0001 13 

\J.Vf\J\J 1 1 -J 


202235_at 


SLC16A1 


solute carrier family 16, member 1 (monocarboxylic acid transporter 1) 


0.000129 


209406_at 


BAG 2 


BCL2-associated athanogene 2 


0.000134 


214577_at 


MAP1B 


microtubule-associated protein 1B 


0.000147 


215404_x_at 


FGFR1 


fibrob ast growth factor receptor 1 


0.000153 


229724_at 


GABRB3 


gamma-aminobutyric acid (GABA) A receptor, beta 3 


0.000183 


215034_s_at 


TM4SF1 


tran^mprnhnnp A I ^iv famik/ mpmhpr 1 
LiQiiJiiicriiiL.'luiitr ^ i_ _h a i g i i i i i y iiicriiikjcrr i 


0 0001 94 


224809_x_at 


TINF2 


TERF1 (TRFI)-interacting nuclear factor 2 


0.000212 


212223_at 


IDS 


iduronate 2-sulfatase 


0.000216 


223144_s_at 


AKIRIN2 


akirin 2 


0.000226 


225267_at 


KPNA4 


karyopherin alpha 4 (importin a pha 3) 


0.000228 


204186_s_at 


PPID 


peptidy prolyl isomerase D 


0.000228 


200895_s_at 


FKBP4 


FK506 binding protein 4, 59 kDa 


0.000269 


226541_at 


FBXO30 


F-box protein 30 


0.000307 


20404 1_at 


MAOB 


monoamine oxidase B 


0.000372 


1 55791 5_s_at 


GST01 


n utathionp S-trrinsfpr^sp ninpns 1 


0 000435 


202199_s_at 


SRPK1 


SFRS nrntpin kinrisp 1 

_ii * \~j ui ULdii i rviiicoc 


0 000476 


211696_x_at 


HBB 


hemoglobin, beta 


0.000491 


224888_at 


SELI 


selenoprotein 1 


0.000593 


227408_s_at 


SNX25 


sorting nexin 25 


0.000618 


222703_s_at 


YRDC 


yrdC domain containing (E. coli) 


0.000654 


229355_at 


NA 


NA 


0.000654 


200065_s_at 


ARF1 


ADP-ribosylation factor 1 


0.000666 


210260_s_at 


TNFAIP8 


tumor necrosis factor, alpha-induced protein 8 


0.000669 


231629_x_at 


KLK3 


kallikrein-related peptidase 3 


0.000683 


201829_at 


NET1 


neuroepithelial cell transforming 1 


0.000719 


240633_at 


DOK7 


docking protein 7 


0.000724 


225766_s_at 


TNP01 


transports 1 


0.00083 


219334_s_at 


OBFG2A 


oligonucleotide/oligosaccharide-binding fold containing 2A 


0.000846 


227622_at 


PCF11 


PCF1 1, cleavage and polyadenylation factor subunit, homolog (S. cerevisiae) 


0.000872 


21 7561_at 


CALGA 


calcitonin-related polypeptide alpha 


0.000891 


214638_s_at 


CCNT2 


cyclin T2 


0.000898 


218624_s_at 


MGC2752 


hypothetical LOG65996 


0.000989 



genes such as G0S2, GSPT1 and BAG5. These activated 
genes represented a stress response for the acute Cd 2+ 
toxicity. 

Our next aim was to identify the gene regulatory path- 
ways that were activated in the HPT cells in response to 
the acute, 1 day exposure to Cd 2+ . The gene pathway 
database compiled in Kyoto Encyclopedia of Genes and 
Genomes (KEGG) has been a valuable source for investi- 
gating pathways. The gene set enrichment analysis 



(GSEA) method was employed to identify activated 
KEGG pathways using the microarray data. The approach 
of GSEA is to test whether the ranks of P values of a set 
of genes are randomly distributed among all the P values 
[12]. Because most of GSEA methods employ two sample 
comparison method such as unpaired or paired t test 
[13], we set out to develop a statistical testing method for 
GSEA based on regression analysis for each individual 
gene. As described earlier, out of a total of 28,720 genes, 
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Table 2 Cd 2+ 


repressed genes in acute cadmium-exposed human proximal tubule cells (1 day, p value < 0.001) 


Probe 


Symbol 


Gene name 


Pr(> |t|) 


213067_at 


MYH10 


myosin, heavy chain 10, non-musce 


4.01 E-05 


225276_at 


GSPT1 


G1 to S phase transition 1 


8.75 E-05 


232349_x_at 


DCAF6 


DHBI and C\ II 4 associated factor 6 


0 000 1 1 5 


1557986_s_at 


SMCR8 


Smith-Magenis syndrome chromosome region, candidate 8 


0.000121 


223141_at 


UCK1 


uridine-cytidine kinase 1 


0.000123 


209196_at 


WDR46 


WD repeat domain 46 


0.000125 


227477_at 


ZMYND19 


zinc finger, MYND-type containing 19 


0.000149 


223396_at 


TMEM60 


transmembrane protein 60 


0.00018 


227765_at 


NA 


NA 


0.000189 


23581 2_at 


TMEM188 


transmembrane protein 188 


0.000237 


239464_at 


NA 


NA 


0.000239 


202337_at 


PMF1 


nnlvflminp-mnrli ilritprl fartor 1 

UUIyul Mil IC II IUUUIOLCU 1 CILLWI 1 


0 000747 


1570506_at 


NA 


NA 


0.000247 


1 555630_a_at 


RAB34 


RAB34, member RAS oncogene family 


0.000249 


22583 1_at 


LUZP1 


eucine zipper protein 1 


0.000254 


204027_s_at 


METTL1 


methyltransferase ike 1 


0.000295 


212694_s_at 


PCCB 


propiony Coenzyme A carboxylase, beta polypeptide 


0.000304 


204104_at 


SNAPC2 


^m^ I niir p^r RMA rirtivstinn mmn py nnlvnpntirlp 7 4-S kD^ 

Jl 1 lull 1 lULIUUI III N/ \ U^LIVOLIIIV^j ^.WIIIL^I^IA, LJ\^J\yLJ'CZ^JI_\\-i\I Z_, \\LyG 


0 000369 


216308_x_at 


GRHPR 


n K/nvvlstp rprli irts^p/hvrlmyvnvn iv^itp rpdiirt^p 


0 000397 


233233_at 


NA 


NA 


0.000438 


200980_s_at 


PDHA1 


pyruvate dehydrogenase (lipoamide) alpha 1 


0.000468 


218042_at 


C0PS4 


C0P9 constitutive photomorphogenic homolog subunit 4 (Arabidopsis) 


0.000476 


210829_s_at 


SSBP2 


single-stranded DNA binding protein 2 


0.000504 


203226_s_at 


TSPAN31 


tetraspanin 31 


0.000516 


202266_at 


TTRAP 


TRAP flnd TNF rprpntnr ^^nri^tpd nmtpin 


0 000541 


221856_s_at 


FAM63A 


family with sequence similarity 63, member A 


0.00066 


227224_at 


RALGPS2 


Ral GEF with PH domain and SH3 binding motif 2 


0000663 


243309_at 


FLJ27352 


hypothetical L0C1 45788 


0.000673 


202883_s_at 


PPP2R1 B 


protein phosphatase 2 (formerly 2A), regulatory subunit A, beta isoform 


0.000684 


201801_s_at 


SLC29A1 


soliitp carrier famiv 79 (nucleoside transnortersl mpmhpr 1 


0 000688 


220974_x_at 


SFXN3 


sideroflexin 3 


0 000699 


1555679_a_at 


RTN4IP1 


reticu on 4 interacting protein 1 


0.000777 


212872_s_at 


MED20 


mediator complex subunit 20 


0.00078 


201888_s_at 


L13RA1 


interleukin 13 receptor, alpha 1 


0.000787 


227567_at 


NA 


NA 


0.00079 


229146_at 


C7orf31 


chromosome 7 open reading frame 31 


0.000806 


224779_s_at 


FAM96A 


family with sequence similarity 96, member A 


0.000843 


21 3022_s_at 


UTRN 


utrophin 


0.000846 


22671 0_at 


C8orf82 


chromosome 8 open reading frame 82 


0.000846 


21 31 81_s_at 


MOCS1 


molybdenum cofactor synthesis 1 


0.000878 


218268_at 


TBC1D15 


TBC1 domain family, member 15 


0.000883 


21 1 708_s_at 


SCD 


stearoyl-CoA desaturase (delta-9-desaturase) 


0.000916 


202299_s_at 


HBXIP 


hepatitis B virus x interacting protein 


0.000948 


213574_s_at 


KPNB1 


karyopherin (importin) beta 1 


0.00097 


205052_at 


AUH 


AU RNA binding protein/enoyl-Goenzyme A hydratase 


0.00097 



4036 were found to have a significant P values (< 0.05) 
correlating to Cd 2+ concentrations. Suppose a KEGG 
pathway has n genes and x genes are significant using 
regression analysis. Our GSEA method consists of deriv- 
ing a hypergemetric probability function shown below. 



4036\ /28720-4036\ 
x J \ n — x J 

22000\ 
n I 
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Using this function, we can test whether x is signifi- 
cantly larger than what would be expected from a 
hypergeometric distribution. Statistically, if we observe r 
genes in the significant set, P(X > r) designates the P 
value for the significance of r genes. Thus, P value = P 
(X > r) = p(r+l) + p(r+2) + p(r+3) + ... +p(x). This gene 
set enrichment analysis found 3 significant pathways, 
NOD-like receptor signaling pathway, base excision 
repair, and steroid biosynthesis, with p values 0.0486, 
0.0385, 0.0009, respectively. NOD-like receptors play an 
important role in regulation of inflammatory and apop- 
totic responses and are activated to stress response [14]. 
DNA repair is a mode of action known to response to 
heavy metal toxicity. Steroid biosynthesis is a metabolic 
process that is required for recovery from an acute toxi- 
city attack. 

Gene expression variation in response to long term 
(13 day) cadmium toxicity 

HPT cells were exposed to Cd 2+ for a long period in order 
to determine the mechanism by which global gene expres- 
sion adapted to environmental challenge. HPT cells were 
exposed to various concentrations of Cd 2+ , 4.5, 9 and 
27 uM Cd +2 for 13 days. Lower concentration of Cd 2+ 
were used for the 13 day time course since 45 uM Cd 2+ 
caused toxicity and a loss of cell viability. There was no 
loss of cell viability when the cells were treated with 4.5 
and 9 uM Cd 2+ over a 13 day period, while exposure to 
27 |iM Cd +2 caused an approximate 30% loss in cell viabi- 
lity [10]. The ability of the HPT cells to form domes was 
only affected by exposure to 27 uM Cd +2 for 13 days. The 
microarray data were processed and analyzed as described 
above for 1 day exposure. The distribution of P values for 
each individual gene is shown in Figure lb. The distribu- 
tion is significantly skewed to small P values, implying a 
large number of genes were altered by cadmium exposure. 
From a total of 28,720 genes, 366 were positively corre- 
lated with the cadmium concentration, and 517 were 
negatively correlated at a P value cutoff of 0.01 (Table 3 
and 4). Similar to that observed for 1 day exposure, long 
term cadmium exposure altered the expression of a num- 
ber of cell cycle control and cell death genes, including 
GSPT1, G2E3 and BCL2. The expression of metal binding 
and transport genes, such as CAB39L, TMC5, SLC5A3 
and SLC35F2, were affected as well. However, expression 
levels of stress response genes such as heat shock and 
DNA repair genes were not changed after 13 days expo- 
sure, instead, many tumor-related genes were activated. 
These included RAB42, DEK, RRAD, TPD52L1 and a vari- 
ety of mitogen related kinases. The GSEA analysis revealed 
9 altered KEGG pathways including 3 cancer pathways, 
prostate cancer, melanoma and p53 signaling pathway 
(Table 5). None of the pathways overlapped with the ones 
identified after 24 hours exposure. 



The common expression patterns between 1 day and 13 
day of cadmium exposure 

In order to investigate genes that showed expression varia- 
tion throughout both short and long term cadmium expo- 
sure, we conducted analysis for the pooled data of the two 
time points. Figure 2 shows a heat map for the 8 samples 
for 1 day and 13 day exposure. The data were normalized 
to a grand mean of 0 and only 901 genes with standard 
deviation > 1 were used for this analysis. Hierarchical clus- 
tering method that was based on Pearson dissimilarity and 
Ward linkage showed as expected that the HPT cells trea- 
ted with high doses of Cd 2+ tended to be grouped 
together. The cells exposed to 27 and 45 uM Cd 2+ for 
both 1 day and 13 days formed a group, and all cells 
exposed to 9 uM or lower Cd 2+ concentration formed the 
other major group. The change of dissimilarity and linkage 
methods did not alter the results. There were 17 genes 
that were positively correlated with cadmium concentra- 
tion for both 1 day and 13 day exposure (Figure 3a). Given 
that the number of positively correlated genes was 403 
and 366 for 1 and 13 day respectively, the random number 
of overlapping genes can be modeled as a hypergeometric 
distribution. The probability of observing more than 17 
overlapping genes from the random sets of 403 and 366 
genes is 5.89E-06. Thus, the two exposure time periods 
had a significant number of common genes that were 
positively correlated with cadmium concentration. Simi- 
larly, we found 26 genes that were negative correlated 
with cadmium concentration for both time periods (Figure 
3b). This number of overlapping genes was statistically sig- 
nificant when compared to selecting 26 gene by chance (P 
value = 1.32E-06). Therefore, we conclude that the HPT 
cells have a common response component in global gene 
expression variation during short term and long term cad- 
mium exposure. 

Gene regulatory network in response to cadmium 
exposure 

In response to heavy metal toxicity, thousands of genes 
will be activated forming a complicated gene regulatory 
network interacting cooperatively. To decode the complete 
gene network for cadmium response is impossible for the 
small sample size used in this study and high complexity 
of network structure [15]. It is also computational forbid- 
den to construct a network by including a large number of 
genes. Therefore, we modeled the gene network focusing 
on only the 901 identified genes as having high variability 
shown in the heat map. Mutual information for transcrip- 
tional levels was used to quantify the relations between 
genes. The resulting adjacent matrix was analyzed for the 
presence of an active network module that has high con- 
nectivity within the module and low connectivity with 
other components of the network. Figure 4 shows a net- 
work modules consisting of 8 genes. DNAJA4 encodes the 
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Table 3 Cd 2+ induced genes in chronic cadmium-exposed human proximal tubule cells (13 days, p value < 0.001) 



Probe 


Symbol 


Gene name 


Pr(> |t|) 


ZZZ>y I 4_S_dl 


f~AR5QI 


r~a \r~\\ i m ta i n *a i ri /a far*at/aiia 5Q i l/ia 

LdiLiuiii uinuiny piuieiii iikc 


1 ^qf n^ 

I .oyt Uj 


910.579 c at 
Z 1 Uj / Z_S_d I 


TPn^9i 1 

I rUjZL I 


tumor protein D52-like 1 


-) 75p r\r: 
Z./jC-Uj 


9 1 1 /1QQ c at 

z i 1 4yy_s_at 


M APl^l 1 
IVIArM I 


mitogen-activated protein kinase 1 1 


5 H7F nq 


->->-> £T Q 7 c nt 

zzzd o/_s_at 


C A I MT7 
bALIN I / 


UDP-N-acetyl-alpha-D galactosamine:polypeptide N-acetylgalactosaminyltransferase 7 (GalNAc-T7) 


D.U3t-UD 


955^90 c at 

Z99DZo_S_at 


uA I jLj 


GATS protein- ike 3 


q q^f nq 

o.ojt-Uj 


91955^ at 
Z I ZjjJ_dl 




glucosamine (N-acety )-6-su fatase 


Q n^ 


9 1 77/I Q -it 

z I / / 4o_at 


ADIrUn 1 


adiponectin receptor 1 


n nnm 1 1 
U.UUU I I z 


999^01 at 


1 IRA£ 


ubipuitin-like modifier activating enzyme 6 


n nnm qq 

U.UUU I jj 


9 1 05/10 c it 

z i y^4y_s_ar 


tAULZ 


exocyst complex component 2 


U.UUU I oo 


91S551Q at 
Z I Oj I y_d L 


pci M 
rLLI I 


y~> a I i la r\ la r\ rva r\ /a /a 1 /Tl r/a c /a la la 1 a ^ 

ptriiinu iiuiiiuKJy i \uiusopi ma,) 


U.UUUZ I 9 


90^^^^ c at 
ZUjjjj i dl_ 


K/lAP4K c ; 
IVIAr 'H-PO 


m itr~iiaian-ar"ti\/atoH tarrvtoin l/iiaaco Uraaco l/iraaco l/inaco ^ 
IIIILUytrM dLUVdLfcJU piULtMII Mlldic Mlldbfc! Mlldbfc! Kllldbc J 


u.uuuz i y 


1 Q£.fM7Q -i -it 

i Dou/ /y_a_at 


M A 

In A 


M A 
IN A 


n nnmi7 
U.UUUzz/ 


9/1/IH^A at 

z44Ujo_at 


jr I AZ 


surfactant associated 2 


n nnn95^ 
U.UUUzdo 


">'~>5707 it 

zzo /y/_at 


rnUzojz 


nypotneticai protein rKUzojz 


n nnm^A 
U.UUUdoo 


9 1 o^on r it 
z i yDou_s_at 


1 IVILD 


transmembrane channel- ike 5 


n nnn'27'2 

U.UUU9/ 9 


1 qq5.567 a at 


LUAODZ 


r\ftf\f~\~\ rnmo r~ laviHaco ci ilai init \/lK raialwiaoiatiHo 9 ('toctic^ 
Ly LULI II Ul lit! L UAlUdit! bUUUIIIL VIU pUlyptipLlUc Z VLcjLlbj 


o nnn.41 ^ 

U.UUUt I D 


9nn^i Q at 
ZUUO I o_at 


1 AQP1 
LAztr I 


1 \kA aia/H nrntoin 1 

liivi diiG jnj proiein i 


n nnn/i 1 /i 

U.UUU4 I 4 


9D5Q9 c ; at 
ZU9yZD_dl 


uL,l_IVI 


glutamate-cysteine igase, modifier subunit 


U.UUU9UZ 


men oq r it 
zuy i yo_s_at 


C^/T1 1 
9 1 I I I 


synaptotagmin XI 


n nnn^no 

u.uuuz>uy 


9151 fl7 at 
Z I j I U/_dl_ 


TMIk 1 
1 INIIV 


i KArz diiu iNL.r\ iriLeidcuny Kinase 


n nnn^ 59 

U.UUU99Z 


9195/11 at 

z i Z34 1 _at 


VIPF^ 
ill rO 


Yip1 domain family, member 6 


n nnn^^f; 

U.UUU909 


9957QQ at 

zzo /yy_at 


i/i A A 1 
IxlAA I oZO 


i/i A A 1 S9f^ 
IxlAA I oZD 


n nnn^QQ 
u.uuujyo 


99QQ1H at 

zzoo i u_at 


rr"MVi 1 


eye in Y-like 1 


n nnn^^f; 

U.UUUOOZ) 


9 1 Q/l 79 c at 
Z I o4/Z_S_at 


ppi n 


pelota homolog (Drosophila) 


n nnn^^Q 
u.uuuooy 


997H.99 at 
ZZ/ UZZ_dL 


rMPHA9 
oIN r UAZ 


glucosamine-6-phosphate deaminase 2 


n nnn^QQ 

U.UUUOoo 


1^9719 a at 
I 00Z/ I Z_d_dl 


MMM AT9 
IMIVIIM A I Z 


ia ir"(ati la a m iHra ni iy"lQ/ati/aja a/aoiawwtraia cfor^cQ 9 

niLOLiNdi ' nue iiucieoLiue dueiiyiyiLidiisieidse z 


u.uuuoy 


9171/1/1 at 
Z 1 / 1 44_d I 


I IRR 
UDD 


ubiguitin B 


n nnn^QQ 
u.uuuoyy 


1559514_at 


L0C1 001 32077 


hypothetical protein L0C1 001 32077 


0.000707 


1 552946_at 


ZNF114 


zinc finger protein 1 14 


0.000794 


22607 1_at 


ADAMTSL4 


ADAMTS-like 4 


0.000822 


204802_at 


RRAD 


Ras-related associated with diabetes 


0.000874 


20881 6_x_at 


ANXA2P2 


annexin A2 pseudogene 2 


0.00088 


212791_at 


C1orf216 


chromosome 1 open reading frame 216 


0.000882 


215704_at 


FLG 


filaggrin 


0.000923 


200638_s_at 


YWHAZ 


tyrosine 3-monooxygenase/tryptophan 5-monooxygenase activation protein, zeta polypeptide 


0.000945 



heat shock protein, SLC7A11 an anion amino acid trans- 
porter, AKR1B10 and AKR1C1 are subunits for a reduc- 
tase associated with cancer, ASPM and KIF20A are 
involved in mitotic spindle functions, ZFAND2A is a zinc- 
finger protein that is inducible by arsenite, and INHBA is 
known to inhibit cell proliferation and to have tumor-sup- 
pressor activity. The network modules provide hypotheses 
for gene interactions which can be verified in subsequent 
biological experiments. 

Discussion 

Microarray experiments have become a standard approach 
to simultaneously quantify the transcriptional levels of 
thousands genes under various conditions. However, due 
to the high cost of chips, microarray experiments are often 



carried out with a limited sample size, which requires a 
careful experimental design. In this study, we chose two 
time points, 1 day and 13 day, to account for toxicity 
caused by short and long cadmium exposure. For each 
time point, we had HPT cells exposed to 3 cadmium con- 
centrations plus a control sample (non-treated). Thus, 
even with only a single replicate at each concentration, we 
were able to test the correlation between transcriptional 
levels and cadmium concentrations. Our analysis suggests 
that experimental design with various conditions over- 
comes the difficulty of lack of replications. 

Exposure of HPT cells to 1 day and 13 days of cad- 
mium induces distinct sets of genes. At 1 day, HPT cells 
appeared to have a stress response that activated a num- 
ber of heat shock and DNA repair genes and pathways. 
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Table 4 Cd 2+ repressed genes in chronic cadmium-exposed human proximal tubule cells (13 days, p value < 0.001) 



Probe 


Svmhnl 


Gene name 


Prf> Itll 
' II' 


9 1 om q c it 

z i ouz9_s_at 


KtU 


peroxisomal D3,D2-enoy-CoA isomerase 


3 ftftt r\Q 


Tm~7Q1 c it 

zuz/o i_s_at 


I MDDCZ 
HNrr9l\ 


inositol polyphosphate-5-phosphatase K 


ft Q 1 c r\Q 
O.O i t-UD 


9 1 <;n^ 3 v n+ 

Z I OUD3_X_at 


P A \fl 1 Q9 A 
rAIVI I ozA 


family with sequence simi arity 182, member A 


n nnn 1 
U.UUU I 


99Qftft^ it 
ZZoDD3_dI 


rvvRi 


cysteine/tyrosine-rich 1 


n nnm no 
u.uuu I uy 


99R374 at 
zzo j / 4_d L 


fi nnrf9R 

V_ I UUI I zo 


LI 1 1 Ul I ILoUl lit! IU Uptrll IfcrdUIIILj lldlllt: ZO 


n 0001 1 q 
u.uuu i i y 


994Q1 3 c at 


i nviivou 


UdllilULdit: Ul lllllfcrl [III LULl lUIIUIIdl IlltrlllUldllt! JU IIUIMUIUy ^J. LtrlLVIildt:^ 


n 0001 3 1 

U.UUU 1 9 1 


9 1 3Q v at 
Z I 3430_X_dI 


ujr I I 


G1 to S phase transition 1 


n nnm 3 9 

U.UUU I 9Z 


9 1 Q/in/i it 
z 1 o4U4_at 


cmvi n 
j IN A I U 


sorting nexin 10 


n nnm 3 3 

U.UUU I 99 


99ftmft at 

zzou i o_at 


LU't/ 


CD47 molecule 


n nnm 3^ 

U.UUU I 99 


9n3 c ;74 at 
ZUjj / 4_dl 


MFII 3 


ni i/~Laar far~t^r iiitarlai il/in 3 rarw \\itar\ 

iiuLiedi idLLOi, iiiLeneuKiu 3 [eyuidLeu 


n 0001 3ft 

U.UUU I 90 


9971 Qn at 
ZZ / I y U_d L 


TMFM37 
I IVILIVI9/ 


transmembrane protein 37 


n nooi 7 

U.UUU I / 


im 11Q e it 

zu I zzD_s_at 


CDDM 1 

onrwl I 


serine/arginine repetitive matrix 1 


n nnn9 1 1 
U.UUUz I I 


I 33Zo4D_S_at 


PAR/19 
nAD4Z 


RAB42, member RAS oncogene family 


n nnn9 1 1 
U.UUUz I / 


99Qc;9n c at 
ZZo3ZU_S_dL 


API P9 
ArLrZ 


diiiyiuiu ueid \r\^t) piecuibui-iiKe pioteiu z 


0 0009 1 Q 

u.uuuz i y 


9977ftQ at 

zz/ /oy_ai 


MA 
InA 


M A 
IMA 


n nnn9 3 1 

U.UUUZ9 I 


ZU 1 oj D_b_dL 


A I IVIIIn 


A 1 IVI 111 Lei dLLOl 


0 0009^ft 
U.UUUZ9D 


91Q37Q v at 

z i y j / y_x_d l 


7MC3CQ 
Z_INi3 jo 


zinc finger protein 358 


0 00097^ 
U.UUUZ / 9 


3 9Q 1 1 it 

ozo i i _at 


IVI Y U I L 


myosin IC 


n nnn9Q9 
U.UUUzoz 


ir\r\~!f\~7 it 


rh\l\Lbn 


protein kinase C substrate 80K-H 


n nnn90ft 
U.UUUzyD 


~)~)C\C\C\1 at 

zzyuu/_dL 


LLJL-Zo3/ OO 


FSHD region gene 1 pseudogene 


n nnn3n 1 

U.UUU9U I 


219271 at 




I I r~)P-l\l-ar"*at\/ -a nhs-fl-n a artricam i n<a-i~ii~i \/ neurit iHia M arph/ na ar~t(~ica minv t ca ncfpracia 1 A. /T^a Kl Ai~-T1 A.\ 
ULT IN dLtrLyl d I pi la U LjdlaL LUidl l 1 1 1 Itr.pUly [Jt:[JLIUt: I NdLtriy ILjd IdLLUid I I III lylLldl liltrl dit: I *+ l,Odl I N nL I \Hj 


0 000304 
u.uuuou^t 


901 Rift c at 
ZU I o I 0_i_d L 


uDAj 


yiiouidsiuiiid diiipinieu seLjueiiLe 


0 00030^ 
U.UUU9U9 


917817 at 
Z I / O I / _d L 


ARPf~4 
Anr LZ-t 


ar"tin ro atQ^ nrntoin 9/3 r~r\ m n La v ciiliiiriit/l 9n 1/ F^a 
dLLIIl leidieu piOieill Z/3 LUllipieX, SUUUIIIL 4, ZU KL/d 


0 0003 1 ^ 
U.UUU9 I 9 


9m 1 Qc; it 
ZU 1 1 OD_ai 


WTD A 1 
n I nA I 


HtrA serine peptidase 1 


n nnn3 1 ft 

U.UUU9 I 0 


9n^730 \/ it 

zU9/3y_X_at 


zJ\lr I U7 


zinc finger protein 107 


n nnn3 93 

U.UUU9Z9 


9 1 9/1 1 Q at 
Z I Z4 I o_a [ 


FI F1 
ELr I 


E74-like factor 1 (ets domain transcription factor) 


n nnn3 9 7 

U.UUU9Z/ 


inQ/1 31 cat 

zuy43 i_s_at 


P AT71 
m I z. I 


POZ (BTB) and AT hook containing zinc finger 1 


n nnn3ft7 

U.UUU9D/ 


9m ^/lO v at 

zu I 34y_x_ai 


ixUIVDD 


ysine (K)-specific demethylase 5B 


n nnn373 

U.UUU9/9 


903R^R <; at 
ZUjOjo b d I 


COX 10 


(^r~lY1H V\rwc\f\\r\c\ r\/tnrhrnmo r nviHaco accomhk/ nrntoin homo A 1 farnocwltrancforaco ('v/oact^ 
LUA 1 U 1 IUI 1 lUIUy , Ly LULl 1 1 Ul 1 It! L UXIUdit! dbitrl IlUly piULtrlll, Mel 1 lc A. Idll Itriy 1 Lldllilfcrldbt! ^ytrdbLj 


0 00049^ 
U.UUUH-Z9 


9 1 QA Q3 c it 

Z I o4o3_s_ai 


r~ 1 1 ^ rf ftn 
L. I I OiTOU 


chromosome 1 1 open reading frame 60 


n c\c\C\aa 
U.UUU44 


9 1 Qn/1 ^ v at 

z i oU4j_x_ai 


r I IVI j 


parathymosin 


n nnn/m/i 

U.UUU494 


9 1 9/I03 r it 

z 1 z4y3_S_at 


bt I Uz 


SET domain containing 2 


U.UUU4jy 


904530^ c at 
ZU4oU 3_S_d L 


1-11 FV 
n I rA 


H1 histone family, member X 


0 00047 1 
U.UUU4/ I 


99ftQ9ft at 

zz oyzo_ai 


UIVlixIN 


dermokine 


n nnn/i i a 

U.UUU4/ 4 


911 mc: v it 

z i i uz9_x_ai 


LUAjd 


cytochrome c oxidase subunit Vb 


n nnn/iQ/i 

U.UUU4y4 


919QnA at 
Z I ZyU4_dl 


I RRrzL7 


eucine rich repeat containing 47 


0 000^44 
U.UUU944 


9 1 ^/l 3 c at 

z i j343_s_ar 


I ARrF 
LA nor 


ike-g ycosy transferase 


n nnn^/ift 

U.UUU94D 


99QQQQ at 

zzoyyo d i 


TMRCfiR 
I IN nv_DD 


tnni \n ontiHo ronoat r^matairnnt/~T ftR 
LIIIIULILULlUt! [LpfcrdL LUI 1 Id III II iy OD 


n 000^7^ 

U.UUU9 / 9 


91 mm at 
z i uzzu_ai 


F7P19 
rz_UZ 


frizzled homo og 2 (Drosophila) 


n nnnft3 9 

U.UUUD9Z 


99^947 at 
ZZ9Z4/_dl 


C 1 On rf ft 

^ i yui ID 


y"Kmrvmc/imo 1 Q /inon roaHiiir~i framo ft 

li ii ui i lusui 1 ie iy upeu leduiuy iidine o 


o oooft^ 

U.UUUD9 


9199AZL c at 
Z I ZZD4_b_d L 


VA/APAI 


wings apart-like homolog (Drosophila) 


0 OOOftft 
U.UUUDD 


9/1/19Q7 it 

z44zo/_at 


CEDC 1 9 
9rK9 I Z 


sp icing factor, arginine/serine-rich 12 


n nnnftft 1 
U.UUUDD I 


991 Q90 r it 

zz 1 ozy_S_ai 


TMDH1 
I INKU I 


transportin 1 


U.UUU/ U4 


1 55331 3_s_at 


SLC5A3 


solute carrier family 5 (sodium/myo-inositol cotrans porter), member 3 


0.000786 


23241 2_at 


FBXL20 


F-box and leucine-rich repeat protein 20 


0.000788 


217964_at 


TTC19 


tetratricopeptide repeat domain 19 


0.000788 


223056_s_at 


XP05 


exportin 5 


0.000792 


213470_s_at 


HNRNPH1 


heterogeneous nuclear ribonucleoprotein H1 (H) 


0.000816 


214543_x_at 


QKI 


quaking homolog, KH domain RNA binding (mouse) 


0.00085 


243829_at 


BRAF 


v-raf murine sarcoma viral oncogene homolog B1 


0.000864 


34260_at 


TEL02 


TEL2, telomere maintenance 2, homolog (S. cerevisiae) 


0.000874 
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Table 4 Cd 2+ repressed genes in chronic cadmium-exposed human proximal tubule cells (13 days, p value ?<? 0.001) 



(Continued) 


Z{JDDyZ_S_dl 


U I Dr I 


C-terminal binding protein 1 


n nnnQ77 
U.UUUo/ / 


OH 1 OTP 3t 




plakophilin 4 


U.UUUooo 


23281 4_x_at 


KLC1 


kinesin light chain 1 


0.000907 


21211 4_at 


LOC552889 


hypothetical protein LOC552889 


0.000922 


225399_at 


TSEN15 


tRNA splicing endonuclease 15 homolog (S. cerevisiae) 


0.000924 


236035_at 


NA 


NA 


0.000931 


229400_at 


HOXD10 


homeobox D10 


0.000961 


208847_s_at 


ADH5 


alcohol dehydrogenase 5 (class III), chi polypeptide 


0.000984 


214800_x_at 


BTF3 


basic transcription factor 3 


0.000992 


20461 9_s_at 


VCAN 


versican 


0.000997 


218826_at 


SLC35F2 


solute carrier family 35, member F2 


0.000998 



This was expected because stress responses that are 
characterized with induced expression of heat shock 
proteins are acute and transient that last for only a few 
hours to a few days [16-20]. This acute response to Cd 2 
+ with the induction of heat shock proteins, such as 
HSP 27, 60, 70 and 90, have been repeatedly verified in 
our previous studies [21-24]. 

The 13 days exposure was used to model the chronic 
response of kidney to environment of cadmium exposure. 
The molecular changes in the HPT cells in response to 
long term cadmium exposure were manifested with activa- 
tion of many oncogenes and cancer pathways. The asso- 
ciation of cancer occurrence and cadmium exposure has 
been established by numerous cohort studies [25-27]. 
Although 13 days of exposure might not exactly represent 
the chronic effect of cadmium on the kidney living envir- 
onment, our experiment suggests that the cell culture sys- 
tem serves as a valid model for studying the long term 
toxicity of cadmium on the kidney. 

Although the two time periods of cadmium exposure 
appeared to induce two distinct toxicity responses, the 
acute and chronic responses, they shared some common 
features in gene expression. The heat map showed that 
the samples were not separated by time of exposure, 
instead, the samples tended to be grouped by doses 

Table 5 Cd 2+ activated KEGG pathways in chronic 
cadmium-exposed human proximal tubule cells (13 days) 



KEGG Pathway P value 



Lysine biosynthesis 


0.000182016 


Sulfur metabolism 


0.001924676 


Glutathione metabolism 


0.001973874 


Caffeine metabolism 


0.004896048 


Protein processing in endoplasmic reticulum 


0.022267231 


Prostate cancer 


0.035237612 


Melanoma 


0.035316953 


Valine, leucine and isoleucine degradation 


0.038038497 


p53 signaling pathway 


0.041271084 



across two time points. All the HPT samples that were 
treated with 27 uM Cd 2+ for both 1 day and 13 days 
formed a cluster. When looking into the genes that were 
correlated with cadmium concentrations in both 1 day 
and 13 days, we found that the number of overlapping 
genes was significantly larger than what was expected 
from random selection. The set of overlapping genes 
included ion transporter proteins, kinases, and transcrip- 
tional factors. These genes should be further investigated 
to reveal their functions in toxicity responses. 

For any biophysical process, thousands of genes inter- 
act to form a network to maintain their biological func- 
tions. There are commonly two statistical approaches to 
study gene interaction, gene set enrichment analysis 
(GSEA) and networks analysis. GSEA has been success- 
fully applied to microarray data to identify the activated 
pathways. However, GSEA has been so far focused on 
comparing different treatment groups, mostly for two 
groups, treatment and control [12,13,28,29]. For this 
study, we developed a GSEA method that can be used for 
ranks of P values from any statistical tests. We have 
showed that our proposed GSEA statistic follows a 
hypergeometric distribution. Thus, the P value of the 
GSEA is available. 

Gene networks analysis is computationally challenging. 
Networks are often highly complicated, consisting of a 
large number of genes that can be fairly dynamic in nat- 
ure. Decoding such a gene network requires researchers 
to quantify the global gene expression under dynamic 
and variable conditions. The current study provides a 
unique opportunity in that microarray data was acquired 
at multiple time points and at various levels of Cd 2+ 
exposure. We used mutual information to account for 
the relations between genes. Mutual information is able 
to measure gene dependency without assuming data dis- 
tributions and linear relationship. However, the mutual 
information-based network models require huge amount 
of computation that is not feasible for modeling all the 
genes. For this reason, we focused on only a set of most 
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Color Key 




-2 0 2 4 

Value 




1d45u.M 
13d 27u-M 
1d27u.M 
1d9uJV1 
13d 9uJV1 
1dO[iM 
13d4.5[iM 
13d OuM 

j.j 

Figure 2 The heat map of global gene expression for all 8 HPT cell samples. Each row represents a HPT cell sample, and each column 
represents a gene. Only the top 901 genes that have high variability across 8 samples are plotted in the heat map. A dendrogram of 
hierarchical clustering is shown on the left side of the map. 



a) 



Repressed Genes 

1 day 13 day 



b) 



Induced Genes 

1 day 13 day 





Figure 3 Venn diagram depicting the number of genes activated in each treatment group. Shown are cadmium correlated genes specific 
for acute (24 h) and chronic (13 days) exposure and genes common to both treatment groups. The number of cadmium correlated genes in 
each group is also broken down into the number of genes induced or repressed by cadmium. 
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AKR1B10 



AKR1C1 



ZFAND2A 




DNAJA4 



SLC7A1 1 



ASPM 



Figure 4 A network module of 8 genes that interact during cadmium toxicity The edge between two genes indicates they are interacting 
directly. 



variable genes across all samples. The network analysis 
revealed gene network modules consisting of a small 
group of genes that intensively interacted with each 
other. The network modules can be verified in molecu- 
lar experiments and they serve as hypothesized models 
for functional relations of genes when cooperating for 
toxicology responses induced by cadmium. 
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